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ABSTRACT 

We use the variational method to obtain approximate analytical expressions for the 
stationary pulselike solutions in birefringent fibers when differences in both phase veloc- 
ities and group velocities between the two components and rapidly oscillating terms are 
taken into account. After checking the validity of the approximation we study how the 
soliton pulse shape depends on its velocity and nonlinear propagation constant. By nu- 
merically solving the propagation equation we have found that most of these stationary 
solutions are stable. 
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1 Introduction 

Pulse propagation in nonlinear birefringent optical fibers is presently an active area of 
research [1-12]. Even "single-mode" fibers support two degenerate modes polarized in 
two orthogonal directions. The degeneracy can be broken introducing deliberately a 
large amount of birefringence through design modifications giving place to the so called 
polarization-preserving fibers. Weak light pulses propagating in these media maintain its 
state of polarization when they are initially polarized along either of the principal axes. 
However at high intensities nonlinear birefringence can change drastically the propaga- 
tion dynamics. Birefringence can also be a residual effect due to imperfections in the 
manufacturing process. 

The equations that govern pulse propagation in nonlinear birefringent pulses have 
been derived by Menyuk [13|. Two main cases of birefringence: i) high and ii) low, have 



been usually considered separately. In the approximation of high birefringence |113| , |I4 
phase and group velocities of both field components are considered to be different. At 
the same time rapidly oscillating nonlinear terms are neglected. On the other hand, the 
approximation of low birefringence takes into account the difference in phase velocities 
between the two linearly polarized components, but neglects their difference in group 
velocities. The full polarization dynamics of solitons in this last approximation has been 



studied in ITlI . 



The general case has been analyzed in just a few papers ^ ^, |TB|. New vectorial 
stationary solitons for this case were first discovered in |Q, and later on the whole family 
of stationary pulselike solutions was numerically determined by Torner et al. f^. These 



families of solitons correspond to mutually trapped pulses that propagate with a common 
group velocity and which therefore can be called walking vector solitons [16|. Stationary 



solutions as singular points of this infinite dimensional dynamical system play a pivotal 
role to understand the propagation dynamics of any arbitrary input pulse. Here we 



generalize the method described in ||Tj| to find accurate analytical approximations of 
the solitary wave solutions for the case of high birefringence. We generalize it for the 
more general case, when no terms are neglected. Our approximate results are compared 
with existing exact ones |jl|, |[ to prove the validity of our approach. We then examine 
the behavior of the two pulse components when the propagation constant (g) or the 
velocity (v) change. We find that the width of the pulses is determined by q whereas 
the amplitudes ratio ai/a2 is determined by v when the energy is higher than a certain 
value. By means of numerical propagation of the stationary solutions found in this way, 
we observe that most of them are stable. 

The rest of the paper is organized as follows. The problem to be addressed by the 
variational method is described in section |^. Our variational approach is developed in 
section |^. In section ^ we present the numerical results and finally in section ^ we briefly 
summarize the main conclusions. 



2 Statement of the problem 

Pulse propagation in a birefringent optical fiber can be described in terms of two non- 
linearly coupled nonlinear Schrodinger equations. In a reference frame traveling along 
the C, axis with the average group velocity, and in normalized units, this set takes the 
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form ||I3|, |15 



= tSUr + -Urr + {\U\^ + A\V\^) U + (3U + BV^U* , 
= -iW, + + (l^P + A\uf) V-PV + BU^V*, 



(1) 



where U and V are the slowly varying envelopes of the two linearly components of the 
field along the x and y axis respectively, ^ is the normalized propagation coordinate, S 
is the inverse group velocity difference, (3 is half the difference between the propagation 
constants, A is the normalized ratio of the nonlinear susceptibilities, B = 1 — A, t is 
the normalized retarded time, and the asterisk denotes complex conjugation. In linearly 
birefringent fibers A = 2/3, as we set here. The set of Eqs.0 has at least three integrals of 
motion 0, i) the action (total energy), ii) the momentum and iii) the Hamiltonian, which 
are a consequence of i)the translational invariance of Eqs.(|^) relative to phase shifts, ii) 
invariance in r and iii) translational invariance in ^. 

Eqs.|l] have two simple linearly polarized pulselike solutions, viz., linearly polarized 
soliton waves along the slow axis 




V = 




and linearly polarized soliton waves along the fast axis 




(3) 



Their corresponding values for the total energy of the pulse are: 




(4) 
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The stability of these solutions has been determined in Ref.[|TT| for 5 = 0, and for arbitrary 
nonzero values of 6 in Ref . IHl . 



Usually Eqs.|l| are written without the last two terms |T3|, as terms with coefficient f3 
can be eliminated by z dependent phase transformations of U and V, and this in turn 
gives the last nonlinear terms rapid phase variations with z which in principle allow their 
neglect. We must however bear in mind that the last terms in Eq.^ are the only ones 
responsible for energy transfer between both polarizations. And although averaging out 
the fast oscillatory terms has proven to be a good approximation to describe most of the 



observed phenomena p, ^, in the picosecond regime, we will go further in retaining 
these terms in the analysis. 



3 Variational approach 



The variational approach when applied to the single Nonlinear Schrodinger Equation 



(NLS) was introduced by Anderson |jT8[. Since then it has been widely used for coupled 
NLS equations. For the specific case of birefringent fibers in the high birefringence ap- 



proximation it has been used for studying dynamical behaviors H, pi 19, 20, 21 1, as well 



as the stationary case [Q. 

We look for stationary solutions moving at a common velocity v in our frame of 
reference. They can be written [|I| as 



[/(r',0 = Pi(r')exp(zgO 



(5) 
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where r' = t — vC, is the common retarded time. Inserting Eqs. ^ into Eqs.[l|, we get a set 
of ordinary differential equations (ODE's) for Pi^2, which reads 

t{6 - v)Pi + {f3- q)Pi + ^Pi + (|Pi|2 + A\P2\^)Pi + BPiPl = 

-i{6 + v)P2 -{(3 + q)P2 + \P2 + (IP2I' + A|Pir)P2 + BPlP; = (6) 

where the overdots indicate derivative respect to r'. Eqs. (P) can be derived (via the Euler- 
Lagrange equations) from the following Lagrangian L = /f^ Cdr', where the lagrangian 
density C is given by 



-t{6 + v) {-P;P2 + P2P;) + 2(g - /3)|Pi 
+2(g + /5)|P2|2 - |Pi|^ - |P2|^ - 2A|PinP2|2 



|2 



-B [Plp; 2 + Pl *p| j . (7) 

We now assume the following ansatz for Pi 2 

Pi = ai sechipi r') exp {i Cir'), i = 1,2. (8) 
Where the variational parameters ai.2, &i,2, and ci^2 are assumed to have real values. The 



above ansatz is inspired from solutions g, and ^, as well as from the work of Pare [|1J . 
We know that the exact solution has a complex phase chirp induced by the four wave 
mixing term (last term in Eqs.|l]), which we neglect for simplicity but also because from 
previous numerical calculations we observed that the phase dependence on r' was mainly 
lineal 0. 

By introducing Eqs.|| into Eqs.0, and applying the Euler-Lagrange equations to the 
averaged Lagrangian one obtains, after some straightforward algebra, a set of six nonlinear 



coupled equations, viz, 

ci = -5 + v + Bbiall{c-) (9) 

C2 = 6 + v-Bb2all{c-) (10) 

■ ^-A^mo{v)-BalM{hM,c^)=0 (11) 



3 bi 3bi b 
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T - ^ + 1^ - ^T^^o(r/) - 5 alb,^M{b,M, C-) = (12) 

T + ^ - - ^?^"^o(r ^) - 5 a? M(6i, 62, c_) = (13) 
3 02 3o2 Oi 

69 2/^2 2an , a? , 1, 9, <9 , , , , 

+ A-^mo r ' - a &2T^M 61, 62, c_ = 14 

3 O2 302 Oi?7 (702 



where 



/oo 
r sin {2xt) sech'^ (bir) sech"^ {b2T)dT (15) 
-00 

c_ = C2 - ci = 25 - ^QBI{c-) (16) 

/oo 
sech'^t sech'^{xt)dt (17) 
-00 

/oo 
cos{2c-t)sech^{bit)sech^{b2t)dt. (18) 
-00 

= g-/5+c?/2+(5-t;)ci, ^2 = g+/3+c^/2-(5+t;)c2, = 61/62, and Q = 2(a2/6i+a2/62) 
In order to simplify the set of equations ^ - |1^, we need to do some assumptions 
and approximations. We first take into account [|l], |^ that the main contribution to c_ 
becomes from 26, while the integral is a minor correction. Therefore, using the first order 



Born approximation [p2] c_ becomes 



c.^26 + ^^QBI{2S) (19) 



Introducing Eq.O into Eqs.O - 14 we get a reduced set of four nonlinear equations. The 



next step must consist in finding accurate analytical expressions for the integrals /, m, 



and M which can not be exactly integrated. Following ref.[|I3| we expand m in a Taylor 



series around 77 = 1 up to second order, and obtain: 



mix] 



m[x] 



4 2, ,2, 7r\, 
3 + 3<^-^' + 3<l-l5)(l 



— X ) 

2 10 



5 - 5(1 + (—)(! 



X] 



(20) 



The same Taylor expansion for M or J gives much more complicated expressions in terms 



of Polygamma functions |]23[, which should be evaluated numerically. Instead of this 
approach we choose to approximate M and / by replacing hi by (61 + &2)/2. In this way 
the integrals can be analytically performed to give: 



M 



167rc_ 
362 



1 + (^^^ 1 csch 



:21) 



dM -lOSvrc? 



365 



■csch 



l-KC^ 



327rc_[62 + 4c2 
365 



-csch 



27rc_ 



327rV [62 + 4c2 
366 



-coth 



27rc_ 



csch 



27rc_ 



(22) 



I{x) 



167rx[62 + a;2 



365 



7!'C0th{ 



2nx, 



4bx 



2x 62 + 4x2 



csch[——\, 



where 6 = 61 + 62. 



(23) 



From here the procedure goes as follow. Given the equation coefficients S, /?, A{= 2/3), 



and B(= 1/3), together with the parameters v and q, we solve the Eqs.pA] -p!4| to obtain 
6i_2, and ai,2- Then ci,2 are obtained from Eqs.|^ |T0|. 
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4 Numerical results 

We have numerically solved the four coupled nonlinear algebraic equations by a Powell 



hybrid method ||2J]. In all cases a unique solution was instantaneoulsy obtained. In 
order to check the validity of the approximations used to solve the integrals we sometimes 
evaluated numerically M, m, and I at each step getting almost identical results but with 
a much higher CPU time consumption. 

Fig.|l| shows the dependence of the soliton energy Q on the nonlinear propagation 
constant q. This diagram was obtained by means of the above described variational 
method. The solid lines represent the solutions given by Eqs. ^. Different sets of almost 
parallel curves correspond to different values of the parameter 5 whilst the value of v 
hardly influences the shape of these curves, fixing solely the minimum allowed value of q. 
These curves coincide totally with those shown in fig. 1(a) of Ref.[^, which were obtained 
exactly by numerically solving the propagation equation for certain values of q whatever 



V was. 



In addition to the solutions approximated by Eq.||, where ai_2 were assumed to be real, 
and therefore represent solutions where U and V are in phase, there exist solutions where 
U and V are 7r/2 out of phase [|I|. They can be obtained directly from our variational 
approach if we do not imposed ai 2 to be real. Alternatively our approach can account 
for them just by changing the ansatz (Eq.^ for Pi). If we choose the following one: 

Pi = aisechipi T')exp{i Cit') 

P2 = a2sech{b2 T')exp{i{c2T' + 7r/2)), (24) 
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and repeat the process, it is easily realized that the only change produced in the nonlinear 
algebraic equations is to replace B by —B. This indicates that the existence of these 
independent solutions is intimately related to the inclusion of the four wave mixing term 
in the propagation equations. 

Figure |^ shows the fraction of energy in the slow mode (Qi/Q) versus the total en- 
ergy (Q) of the two-parameter family of solutions for 5 = /? = 1 as obtained from our 
variational method. In (a) and (b) Pi is given by Eq.^ and Eq. ^ respectively. The results 
can be compared with the exact ones in Fig. 1 of Ref . [|l| , which were found numerically. 
The symmetry of these curves is a result of the symmetry that possesses the Lagrangian 
under: 

q *q = q^ , V ,v = — -v, 

2 

2(3 

C2,i — > C2,i = — - ci,2, 02 < — > ai, 62 < — ^ hi (25) 


Figs. |T], and ^ illustrate the accuracy of the variational approach as developed in the 
section ^. We must remark that solving the set of four nonlinear algebraic equations is an 
easy task that can be done really fast. In our computer (Alpha Dec 2100/500) we obtain 
thousands of these variational solutions in just a few CPU seconds. 

Fig. shows the variation of the width ratio [t] = 61/62) vs. the propagation constant 
q for the stationary solutions when U and V are (a) in phase, and (b) in quadrature. 
We take (3 = 6 = 1, and three values of v are considered, which are written close to the 
corresponding curves. In all cases rj tends to 1 as g increases, i.e. for high values of Q 
(see Fig. |1]). The departure of t] from unity is higher when U and V are in quadrature 
than when they are in phase. In any case it is small except around the minimum allowed 
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value of g, where the curves emerge and one of the components {V in the cases of figure]^) 
is almost zero. This figure serves also to verify a posteriori that our approximations of 
integrals I, m and M had good basis. 

Similarly Fig.§ shows (a) C2 — Ci and (b) C2 vs. q for the same values of the parameters. 
As expected, at the minimum value allowed for q, where the solution is a fast soliton 
(Eq.|^), C2 is exactly 6 + v. As it occurred for Q vs. q, the value of c_ hardly depends 
on V, being mainly determined by the value of q. It is also remarkable that while the 
solutions for U and V in phase decrease its frequency difference (i.e. c_) as we increase 
q (and therefore Q), the opposite happens for the solutions with U and V in quadrature. 
In these last cases c_ becomes more different to 26 than when U and V are in phase. The 
same happened with 77 respect to 1 (see Fig. ^. In fact we observe that the variational 
solutions were much more accurate for the solutions with U and V in phase than for those 
in quadrature. 

Fig.H shows the variation of the widths and peak amplitudes vs. q for the solutions 
with U and V in phase. It can be seen that for a given value of 6 and (3, the widths 
depend almost exclusively on q, being almost completely independent on v. On the other 
hand the asymptotic behavior of 01/02 does not depend on q but only on v (see Fig.^). 

We have numerically solved Eqs.|l], taking as initial conditions the variational solutions 
corresponding to different values of the parameters v and q and fixed values of 5 and 
/?, viz. (3=1 = 6. In all cases that we propagate in-phase solutions, we obtained 
stable propagation, whilst 7r/2-dephased solutions were sometimes unstable. Fig. ^ shows 
the stable propagation of the in-phase variational solution for q = 3.2 and v = 0.9, 
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whereas Fig.^ shows the same for the 7r/2-dephased solution corresponding to g = 6 and 
V = 0.9. Contrary to the previous case, a small difference can be appreciated between the 
variational solution and the stationary one. The exact stable stationary solution (which 
of course is 7r/2-dephased) is reached after a very short propagation distance. Finally 
Fig. H illustrates the unstable behavior of a 7r/2-dephased solution (g = 7.2, v = 0.8). 
The solution oscillates around the stationary solution with increasing amplitude emitting a 
small quantity of radiation. Eventually, after the emisssion of a large quantity of radiation, 
a stable solution is reached (not shown in the figure). In Figs. r" = r + 5. 



5 Summary 

We have developed an accurate variational approach to derive analytical approximations 
of the coupled pulselike stationary solutions in birefringent fibers in the most general 
case when differences in both phase velocities and group velocities between the two com- 
ponents and fastly oscillatory terms are taken into account (they are important in the 
subpicosecond regime). As particular cases it includes those where any term in Eqs.|I] can 
be neglected. We have shown that the difference between the central frequencies of the 
components, their widht and their energies are almost independent on v, being mainly 
determined by q. We have also shown that the ratio ai/a2 is mainly determined by v. 
The stability of these solutions has been briefly considered, by numerically propagating 
them. In all the cases we propagated solution in-phase, they happened to be stable whilst 
those 7r/2 out of phase present some intervals of stability. A global stability analysis of 
these solutions remains to be done. This variational method could be useful to design 
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soliton-dragging logic gates ^5 . 
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Figure 1: Dependence of the soliton energy Q on the nonhner propagation constant q for 
the fast and slow hnearly polarized solitons (continuous lines) and coupled soliton states 
for /9 = 1, 5 = 2 and 4, and several values of the parameter which is written close to 
the point where the curves emerge. This point is marked with a filled circle 
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Figure 2: Fraction of energy carried by the walking soliton in the slow polarization vs. 
the total energy of the solution, (a) U and V are in phase, (b) U and V have a relative 
phase difference of 7r/2. They correspond to 5 = /3 = 1, and different values of which 
is written in the figure. 




Figure 3: Width ratio (rj) vs q for S — P — 1, and three values of the parameter v. a) U 
and V are in phase, b) U and V are in quadrature. 
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Figure 4: a) Central frequency difference between the soliton components V and U vs. 
q. b) Frequency shift of the fast component as a function of the propagation constant q. 
The dashed hues represent the solutions where U and V are in phase and the dotted ones 
those in quadrature. The values of the parameters are the same than in Fig. ^. 
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Figure 5: a) Inverse widths (&i,2), b) peak amplitudes (01,2) and c) amplitude ratio 
(01/02) of the variational solution for U and V in phase vs. q. In (a) and (b) the dotted 
hnes are for U and the dotted lines for y. S — P — 1 and the values of v are written in 
the figures 
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Figure 6: Stable propagation of a walking soliton with its components U and V in phase. 
p = S = l,q = 3.2, and v = 0.9 
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Figure 7: Stable propagation of a walking soliton with its components U and V in 
quadrature. /? = 5 = 1, g = 6, and v = 0.9 
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Figure 8: Unstable propagation of a stationary solution with its components in quadra- 
ture. (3 = 5 = 1, q = 7.2, and v = 0.8 
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